#!/usr/bin/perl
# 2014/07/03 Co-localization test of specific atribute

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;


use threads;
use Thread::Queue;
use Thread::Semaphore;
use threads::shared;

if(@ARGV != 6 or $ARGV[0] eq '--help'){
	die "Usage : $0 -a [attribute file] -N [total association number] -o [simulation output]\n";
}

my %opt;
getopts("a:N:o:", \%opt);
my $FILE_attribute = $opt{a};
my $TOTAL_ASSOCIATION_NUMBER = $opt{N};
my $FILE_output : shared;
$FILE_output = $opt{o};
#MACS_peak_1
#MACS_peak_2	conv_class1	SPNCRNA	LTR	conv_other
#MACS_peak_3	high	conv_class1	SPNCRNA	LTR	conv_other
#MACS_peak_4	conv_class1	SPNCRNA	conv_other	5SrRNA
#MACS_peak_5	high	conv_class2	conv_other
#MACS_peak_6	SPNCRNA	conv_class2	SPSNORNA	conv_other
#MACS_peak_7	conv_class1	LTR	conv_other


my $NUMBER_OF_SIMULATIONS : shared;
$NUMBER_OF_SIMULATIONS = 10000;
my $THREAD_NUMBER = 50;


#---------------------------------------
# read attribute
#---------------------------------------
my %atr : shared;
my %atr_variation;
my @atrs : shared;
my @binding_sites : shared;
my $fh_at = IO::File->new($FILE_attribute) or die "cannot open $FILE_attribute: $!";
while($_ = $fh_at->getline()){
	s/\r?\n//;
	my ($name, @other) = split /\t/;
	if(@other > 0){
		$atr{$name} = join("\t", @other);
		foreach my $at(@other){
			$atr_variation{$at} = 1;
		}
	}
	push @binding_sites, $name;
}
$fh_at->close();
foreach my $key(keys %atr_variation){
	push @atrs, $key;
}

#---------------------------------------
# set simulation number
#---------------------------------------
my $semaphore_SIM_NUM = Thread::Semaphore->new(1);
my $semaphore_ON_OFF = Thread::Semaphore->new(1);

sub simulationNum{
	if($NUMBER_OF_SIMULATIONS > 0){
		$semaphore_SIM_NUM->down();
		$NUMBER_OF_SIMULATIONS -= 1;
		$semaphore_SIM_NUM->up();
		return 1;
	}else{
		return 0;
	}
}



#---------------------------------------
# output title
#---------------------------------------
my @title : shared;
for(my $i = 0; $i < @atrs; $i++){
	my $at1 = $atrs[$i];
	for(my $j = $i; $j < @atrs; $j++){
		my $at2 = $atrs[$j];
		my $t = $at1 . ":" . $at2;
		push @title, $t;
	}
}
my $fh_out = IO::File->new($FILE_output, 'w') or die "cannot write $FILE_output: $!";
$fh_out->print(join("\t", @title) . "\n");
$fh_out->close();

#---------------------------------------
# simulation
#---------------------------------------
my @threads;
for(1..$THREAD_NUMBER){
	my $thread = new threads(\&simulation);
	push @threads, $thread;
}

sub simulation{
	while(&simulationNum()){

		### make random pairs
		my $TOTAL_PAIRS = 0;
		my %random_pairs;
		while($TOTAL_PAIRS < $TOTAL_ASSOCIATION_NUMBER){
			my $rn1 = int(rand(@binding_sites));
			my $rn2 = int(rand(@binding_sites));

			if($rn1 == $rn2 or exists $random_pairs{"$rn1\t$rn2"} or exists $random_pairs{"$rn2\t$rn1"}){
				next;
			}else{
				$random_pairs{"$rn1\t$rn2"} = 1;
				$TOTAL_PAIRS++;
			}
		}

		### count combination of properties
		my %count_properties;
		foreach my $key(keys %random_pairs){
			my ($rn1, $rn2) = split /\t/, $key;
			my @at1;
			my @at2;
			if(exists $atr{$binding_sites[$rn1]}){
				@at1 = split /\t/, $atr{$binding_sites[$rn1]};
			}
			if(exists $atr{$binding_sites[$rn2]}){
				@at2 = split /\t/, $atr{$binding_sites[$rn2]};
			}
			if(@at1 > 0 and @at2 > 0){
				foreach my $k1(@at1){
					foreach my $k2(@at2){
						$count_properties{"$k1\t$k2"}++;
					}
				}
			}
		}

		my @output;
		foreach my $key(@title){
			my ($at1, $at2) = split /:/, $key;
			my $num = 0;
			if(exists $count_properties{"$at1\t$at2"}){
				$num += $count_properties{"$at1\t$at2"};
			}
			if(exists $count_properties{"$at2\t$at1"}){
				$num += $count_properties{"$at2\t$at1"};
			}
			push @output, $num;
		}

		$semaphore_ON_OFF->down();
		open OUT, ">>$FILE_output" or die "cannot open $FILE_output: $!";
		print OUT join("\t", @output) . "\n";	# output result
		close OUT;
		$semaphore_ON_OFF->up();
	}
}

foreach (@threads){
	$_->join();
}


